In [1]:
# This changes the current directory to the base saga directory - make sure to run this first!
# This is necessary to be able to import the py files and use the right directories,
# while keeping all the notebooks in their own directory.
import os
import sys
if 'saga_base_dir' not in locals():
saga_base_dir = os.path.abspath('..')
if saga_base_dir not in sys.path:
os.chdir(saga_base_dir)
In [2]:
from __future__ import print_function, division
import hosts
import numpy as np
import healpy
from astropy import units as u
from astropy.coordinates import SkyCoord, UnitSphericalRepresentation, CartesianRepresentation
In [3]:
%matplotlib inline
from matplotlib import pyplot as plt
plt.rcParams['image.cmap'] = 'viridis'
In [4]:
hscmap, hschead = healpy.read_map('S16A_wide2_fdfc_i_limitmag_hp_frac_cut_0.02.fits', h=True)
hschead = dict(hschead) #This is a VERY foolish healpy design choice!!
In [5]:
healpy.mollview(hscmap)
In [6]:
def query_disc_astropy(cen, rad, nside=hschead['NSIDE'], **kwargs):
uvec = cen.represent_as(UnitSphericalRepresentation).represent_as(CartesianRepresentation).xyz
return healpy.query_disc(nside, uvec, rad.to(u.radian).value)
In [7]:
scs = SkyCoord([0, 0, 45, 270]*u.deg, [0, 12, 0,45]*u.deg)
rads = [5, 10, 30, 45]*u.deg
pixin = np.zeros_like(hscmap)
for i, (sc, rad) in enumerate(zip(scs, rads)):
pixin[query_disc_astropy(sc, rad)] += i+1
healpy.mollview(pixin)
In [8]:
hostlst = hosts.get_saga_hosts_from_google() #'named' hosts
hostdct = {host.name:host for host in hostlst}
In [10]:
nametofrac_environs = {}
nametofrac_1deg = {}
for host in hostlst:
cen = host.coords
rad = host.environsarcmin*u.arcmin
host_pix = query_disc_astropy(cen, rad)
nametofrac_environs[host.name]= np.sum(hscmap[host_pix])/len(host_pix)
host_pix = query_disc_astropy(cen, 1*u.deg)
nametofrac_1deg[host.name]= np.sum(hscmap[host_pix])/len(host_pix)
nametofrac_environs, nametofrac_1deg
Out[10]:
In [11]:
for hnm, frac in nametofrac_environs.items():
if frac>0:
print(hnm)
print(frac)
print(hostdct[hnm].coords)
print('http://legacysurvey.org/viewer-dev?zoom=10&ra={0.ra.value}&dec={0.dec.value}'.format(hostdct[hnm].coords))
print('')
In [71]:
healpy.cartview(hscmap)
ras = np.array([h.coords.ra.wrap_at(180*u.deg).degree for h in hostlst])
decs = np.array([h.coords.dec.degree for h in hostlst])
overlap = np.array([nametofrac_environs[h.name]>0 for h in hostlst])
plt.gcf().axes[0].scatter(ras[overlap], decs[overlap], c='w', lw=0, s=30)
plt.gcf().axes[0].scatter(ras[~overlap], decs[~overlap], c='k', lw=0, s=30);
In [74]:
healpy.cartview(hscmap)
ras = np.array([h.coords.ra.wrap_at(180*u.deg).degree for h in hostlst])
decs = np.array([h.coords.dec.degree for h in hostlst])
overlap = np.array([nametofrac_environs[h.name]>0 for h in hostlst])
plt.gcf().axes[0].scatter(ras[overlap], decs[overlap], c='w', lw=0, s=50)
plt.gcf().axes[0].scatter(ras[~overlap], decs[~overlap], c='k', lw=0, s=50)
plt.xlim(-180,-90)
plt.ylim(-20,20);
In [81]:
for h in hostlst:
if nametofrac_environs[h.name]>0:
print(h.coords.to_string('hmsdms', sep=':'))
Visual inspection with HSCmap reveals... only Alice actually overlaps! These must be the tracts, not the actual imaging